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ABSTRACT 


In  recent  years,  there  has  been  a  resurgence  of  interest  in  US  Air  Force  in  space  access  and 
therefore  Air  Force  is  interested  in  the  hypersonic  aerodynamics  of  its  future  space  operations 
vehicles,  long-range-strike  vehicles  and  military-reusable  launch  vehicles.  Hypersonic  flows 
about  space  vehicles  produce  flow  fields  in  thermodynamic  non-equilibrium  with  local  Knudscn 
numbers  Kn  -kt L  (where  k  is  the  mean  free  path  of  gas  molecules  and  L  is  a  characteristic 

length)  which  may  lie  in  all  the  three  regimes  -  continuum,  transition  and  rarefied.  Flows  in 
continuum  regime  can  be  modeled  accurately  by  the  Navier-Stokes  (NS)  equations;  however  the 
flows  in  transition  and  rarefied  regimes  require  a  kinetic  approach  such  as  the  Direct  Simulation 
Monte  Carlo  (DSMC)  method  or  the  solution  of  the  Boltzmann  equation. 

The  objective  of  this  research  project  has  been  to  develop  a  computational  methodology  and 
a  code  for  computing  hypersonic  non-equilibrium  shock  wave  flows  of  multi-component  reactive 
gas  mixtures  of  diatomic  gases  using  the  Generalized  Boltzmann  Equation  (same  as  the  Wang- 
Chang  Uhlenbcck  equation  which  accounts  for  the  degenerate  energy  levels)  at  Knudscn 
numbers  in  transitional  and  rarefied  flow  regimes.  Several  milestones  have  been  achieved: 

1.  A  3D  Generalized  Boltzmann  Equation  (GBE)  solver  has  been  developed  for  a  Cartesian 
mesh.  The  solver  has  been  validated  by  computing  the  ID  shock  structure  in  nitrogen  for 
Rotational-Translational  (R-T)  relaxations  and  comparing  the  numerical  results  with  the 
experimental  data  for  Mach  numbers  up  to  15.  The  solver  has  been  exercised  successfully  for 
computing  the  2D  blunt  body  flows  in  nitrogen  and  3D  flow  from  a  rectangular  jet  of  nitrogen  in 
vacuum  for  RT  relaxations.  The  issues  of  stability  of  the  algorithm  and  the  possibility  of 
reducing  the  number  of  rotational  levels  in  the  computations  without  compromising  the  accuracy 
of  the  solutions  have  been  rigorously  addressed. 

2.  A  computational  methodology  has  been  developed  to  compute  the  hypersonic  shock 
structure  in  diatomic  gases  including  both  the  RT  and  Vibrational-Translational  (V-T) 
relaxations.  1-D  shock  structure  in  nitrogen  has  been  computed  including  both  R-T  and  V-T 
relaxations  and  has  been  validated  by  comparing  the  results  with  the  experimental  data. 

3.  A  computational  methodology  has  been  developed  to  compute  the  hypersonic  shock 
structure  in  a  non-reaetive  mixture  of  two  diatomic  gases.  1-D  shock  structure  has  been 
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computed  in  an  inert  mixture  of  nitrogen  and  oxygen  for  R-T  relaxations.  To  accomplish  this,  the 
GBE  is  formulated  and  solved  in  “impulse  space"  instead  of  velocity  space. 

4.  A  new  two-level  kinetic  model  has  been  developed  for  computing  the  RT  relaxations  in  a 
diatomic  gas  and  has  been  validated  by  comparing  the  results  with  the  solutions  of  complete 
GBE.  The  model  is  about  twenty  times  more  efficient  than  the  GBE  in  computing  the  shock 
structure.  It  should  be  noted  that  the  model  is  different  than  the  BGK  model;  it  accounts  for  both 
elastic  and  inelastic  collisions. 
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1.  INTRODUCTION 


In  recent  years,  there  has  been  a  resurgence  of  interest  in  US  Air  Force  in  space  access  and 
therefore  Air  Force  is  interested  in  the  hypersonic  aerodynamics  of  its  future  space  operations 
vehicles,  long-range-strike  vehicles  and  military-reusable  launch  vehicles.  Hypersonic  flows 
about  space  vehicles  produce  flow  fields  in  thermodynamic  non-equilibrium  with  local  Knudsen 
numbers  Kn  -AIL  (where  A  is  the  mean  free  path  of  gas  molecules  and  L  is  a  characteristic 

length)  which  may  lie  in  all  the  three  regimes-  continuum,  transition  and  rarefied.  Therefore, 
there  is  an  important  need  for  a  single  unified  Computational  Fluid  Dynamics  (CFD)  code  that 
can  treat  all  the  three  flow  regimes  in  thermodynamic  non-equilibrium  accurately  and  efficiently. 
Flows  in  continuum  regime  can  be  modeled  accurately  by  the  Nav  ier-Stokes  (NS)  equations; 
however  the  flows  in  transition  and  rarefied  regimes  require  a  kinetic  approach  such  as  the  Direct 
Simulation  Monte  Carlo  (DSMC)  method  or  the  solution  of  the  Boltzmann  equation. 

One  of  the  critical  issues  in  accurate  prediction  of  non-equilibrium  flows  is  the  ability  to 
simulate  the  translational  and  internal  energy  mode  relaxation  of  polyatomic  (in  particular 
diatomic)  molecules  present  in  these  flows.  Relaxation  of  diatomic  molecules  in  non-equilibrium 
flows  is  very  different  from  that  of  monoatomic  molecules  due  to  the  internal  degrees  of 
freedom;  therefore  it  is  important  to  study  the  effect  of  the  internal  degrees  of  freedom  upon  the 
energy  transfer  between  colliding  diatomic  molecules.  It  turns  out  that  the  simulation  of  internal 
energy  mode  relaxation  is  fundamentally  different  in  the  continuum  (NS)  and  kinetic  approaches. 
In  the  continuum  approach,  NS  equations  contain  the  source  terms  of  reaction  probabilities  for 
quantifying  the  thermal  and  chemical  non-equilibrium  effects  which  are  typically  available  from 
experiments  for  equilibrium  conditions  that  have  the  translational  temperature  dependence.  For 
flows  with  Kn  ~  0.01,  this  approach  based  on  NS  equations  is  very  effective  in  computing 
hypersonic  flows  with  small  deviation  from  translational  non-equilibrium  [1],  However,  at 
higher  Kn  for  flows  in  transition  and  rarefied  regimes,  the  kinetic  methods  based  on  the 
Boltzmann  equation  provide  more  detailed  information  on  the  degree  of  non-equilibrium. 

During  the  past  fifteen  years,  following  Bird  [2],  DSMC  methods  have  been  developed  for 
computing  non-equilibrium  flows  of  monoatomic  and  diatomic  gases  [3-5],  Typically,  in  most  of 
the  DSMC  solvers,  the  diatomic  molecules  are  modeled  assuming  quantized  rigid  rotors  for 
rotational  energy  levels  and  anharmonie  oscillators  for  vibrational  energy  levels.  Elastic  eross- 
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sections  are  based  on  Variable  Hard  Sphere  (VHS)  model  and  inelastic  eross-seetions 
(Rotational-Translational  (R-T)  and  Vibrational-Translational  (V-T))  are  based  on  Borgnakke- 
Larsen  model  [6]  assuming  constant  or  temperature  dependent  collision  numbers  ZR. 
Dissociation  cross-sections  are  based  on  the  Weak  Vibrational  Bias  model  [7]  or  its  variants. 
However,  it  has  been  shown  that  the  non-equilibrium  rarefied  flows  of  diatomic  gases,  in  which 
the  gas  moleeules  transfer  energy  among  translational,  rotational  and  vibrational  degrees  of 
freedom,  cannot  be  aeeurately  predicted  by  using  the  simple  collision  models  in  DSMC  methods 
[8].  In  recent  years,  considerable  effort  has  also  been  devoted  toward  the  development  of 
approaches  using  simplified  models  of  Boltzmann  equation  (e.g.  BGK  type  models)  that  include 
the  multi-translational  temperatures,  rotational  relaxation,  and  dissociation  kinetics  [9,10],  which 
have  shown  some  promise.  However,  these  approaches  also  have  limitations,  especially  in 
prediction  of  strong  shocks  encountered  in  hypersonic  non-equilibrium  flow's.  Additionally,  the 
transition  and  rarefied  regimes  are  characterized  by  the  formation  of  narrow  boundary  layers 
with  sharp  variation  in  flow  parameters,  the  zones  with  considerable  compression  of  a  gas  at  the 
seale  of  the  molecular  mean  free  path,  and  the  low  density  stagnation  zones.  Thus  the  most 
accurate  description  of  the  physics  of  these  flows  can  be  obtained  by  solving  the  Boltzmann 
equation  for  a  diatomic  gas,  namely  the  Wang-Chang-Uhlenbeek  Master  Equation  [11]  or  a 
generalized  Boltzmann  equation  for  a  reactive  mixture  of  gases. 

In  solving  the  Boltzmann  equation  by  a  finite-difference  method,  the  principal  difficulties 
arise  in  calculation  of  the  multi-dimensional  collision  integral;  the  approximation  to  the  collision 
integral  must  tend  to  the  actual  one  as  the  mesh  size  in  the  velocity  space  tends  to  zero.  In  recent 
years,  there  has  been  significant  progress  toward  the  development  of  an  efficient  and  accurate 
numerical  method  for  the  solution  of  Boltzmann  equation  for  a  monoatomie  gas  [12-15].  In  this 
method,  the  Boltzmann  equation  is  solved  on  fixed  space  and  velocity  grids  by  a  finite-difference 
method.  A  projection  method  (that  ensures  that  the  velocities  before  and  after  collision  belong  to 
the  same  grid  of  discrete  ordinates)  is  employed  for  the  evaluation  of  the  collision  integral  that 
ensures  exact  conservation  laws  for  mass,  momentum,  and  energy  as  well  as  zero  value  of  the 
integral  under  thermodynamic  equilibrium  (when  the  distribution  function  is  Maxwellian).  The 
last  property  eliminates  the  numerical  error  of  computing  the  principal  part  of  the  solution 
outside  the  Kundsen  layers  and  shock  waves  and  thus  considerably  increases  the  accuracy  and 
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efficiency  of  the  method.  The  differential  part  of  the  Boltzmann  equation  is  approximated  by  an 
explicit  second-order  flux -conservative  scheme.  The  combined  system  of  difference  equations 
(for  the  collision  integral  and  the  differential  part)  is  solved  by  the  splitting  method  which  splits 
the  solution  process  in  two  stages:  the  collision  relaxation  and  free  molecular  flow.  This  method 
has  been  developed  by  Professor  Chcrcmisin  of  the  Computing  Center  of  the  Russian  Academy 
of  Science  [12-15].  It  has  been  extensively  applied  by  him  and  many  other  researchers  [16] 
including  the  author  of  this  report  [17-21],  The  key  numerical  features  of  this  method  arc:  (a)  it  is 
fully  conservative,  (b)  it  preserves  the  positivcncss  of  the  solution,  (c)  it  docs  not  disturb  the 
thermodynamic  equilibrium  and  therefore  can  be  applied  for  computing  flows  approaching 
continuum  regime,  (d)  it  is  essentially  deterministic  and  therefore  docs  not  produce  statistical 
noise,  (c)  it  employs  numerically  efficient  integration  grids  that  make  it  very  efficient,  (0  it  can 
employ  a  variable  mesh  that  may  exceed  the  local  mean  free  path  in  the  regions  of  low  gradients, 
and  (g)  the  method  can  be  easily  parallelized.  It  has  several  advantages  over  the  DSMC  method 
c.g.  the  DSMC  method  requires  mesh  spacing  less  the  mean  free  path  in  the  entire  field,  it 
employs  not  very  realistic  molecular  potentials  (like  VHS)  instead  of  the  more  accepted  ones  like 
the  Lcnnard-Joncs  potential  with  established  parameters  for  each  gas  (c.g.  N2  and  02),  and  for 
inelastic  collisions  DSMC  method  employs  models  that  arc  not  physically  justifiable.  For 
example  in  the  most  commonly  used  Borgnakkc-Larson  model  [6]  in  the  DSMC  method,  the 
molecules  are  divided  in  two  parts:  the  molecules  in  major  part  collide  elastically  and  in  the  rest 
with  internal  -translational  energy  transfer  that  presumes  a  thermodynamic  equilibrium.  This 
model  is  therefore  not  very  accurate. 

The  objective  of  this  research  project  has  been  to  develop  a  computational  methodology  and 
a  code  for  computing  hypersonic  non-equilibrium  shock  wave  flows  of  multi-component  reactive 
gas  mixtures  of  diatomic  gases  using  the  Generalized  Boltzmann  Equation  (same  as  the  Wang- 
Chang  Uhlcnbeck  equation  which  accounts  for  the  degenerate  energy  levels)  at  Knudsen 
numbers  in  transitional  and  rarefied  flow  regimes.  It  should  be  noted  that  in  the  GBE,  the  internal 
and  translational  degrees  of  freedom  arc  considered  in  the  framework  of  quantum  and  classical 
mechanics  respectively.  The  general  computational  methodology  for  the  solution  of  the  GBE  is 
similar  to  that  for  the  classical  BE  for  a  monoatomic  gas  except  that  the  evaluation  of  the 
collision  integral  becomes  significantly  more  complex  due  to  the  quantization  of  rotational  and 
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vibrational  energy  levels.  The  transition  probabilities,  elastic  and  inelastie  eross-seetions  etc.  of  a 
gas  molecule  are  needed  for  the  solution  of  the  collision  integral.  Lennard  -Jones  potential  with 
two  free  parameters  is  used  to  obtain  the  elastic  eross-seetion  of  the  gas  molecules,  and  the  so 
called  “combinatory  relations’'  are  used  to  obtain  parameters  of  Lennard-Jones  potential  for  an 
interaction  of  molecule  A  with  molecule  B  knowing  the  parameters  of  A  and  B  [11].  The 
probability  of  transition  in  inelastie  collisions  is  determined  using  the  approach  by  Beylieh  [22, 
23].  These  inputs  allow  for  the  calculation  of  the  Boltzmann  Collision  Integral  in  GBE  for  a 
diatomic  gas  and  a  reactive  mixture  of  gases. 

Several  milestones  have  been  achieved: 

1.  A  3D  Generalized  Boltzmann  Equation  (GBE)  solver  has  been  developed  for  a  Cartesian 
mesh.  The  solver  has  been  validated  by  computing  the  ID  shock  structure  in  nitrogen  for 
Rotational-Translational  (R-T)  relaxations  and  comparing  the  numerical  results  with  the 
experimental  data  for  Mach  numbers  up  to  15.  The  solver  has  been  exercised  successfully  for 
computing  the  2D  blunt  body  flows  in  nitrogen  and  3D  flow  from  a  rectangular  jet  of  nitrogen  in 
vaeuum  for  RT  relaxations.  The  issues  of  stability  of  the  algorithm  and  the  possibility  of 
reducing  the  number  of  rotational  levels  in  the  computations  without  compromising  the  accuracy 
of  the  solutions  have  been  rigorously  addressed. 

2.  A  computational  methodology  has  been  developed  to  compute  the  hypersonic  shock 
structure  in  diatomic  gases  including  both  the  RT  and  Vibrational-Translational  (V-T) 
relaxations.  1-D  shock  structure  in  nitrogen  has  been  computed  including  both  R-T  and  V-T 
relaxations  and  has  been  validated  by  comparing  the  results  with  experimental  data. 

3.  A  computational  methodology  has  been  developed  to  compute  the  hypersonic  shock 
structure  in  a  non-reaetive  mixture  of  two  diatomic  gases.  1-D  shock  structure  has  been 
computed  in  an  inert  mixture  of  nitrogen  and  oxygen  for  R-T  relaxations.  To  accomplish  this,  the 
GBE  is  formulated  and  solved  in  “impulse  space”  instead  of  velocity  space. 

4.  A  new  two-level  kinetic  model  has  been  developed  for  computing  the  RT  relaxations  in  a 
diatomic  gas  and  has  been  validated  by  comparing  the  results  with  the  solutions  of  complete 
GBE.  The  model  is  about  twenty  times  more  efficient  than  the  GBE  in  computing  the  shock 
structure.  It  should  be  noted  that  the  model  is  different  than  the  BGK  model;  it  accounts  for  both 
elastic  and  inelastie  collisions. 
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2.  TECHNICAL  APPROACH 


As  mentioned  before  in  Section  1,  gas  flows  in  continuum-transition  regime,  where  the 
Knudsen  number  Kn  is  0(1)  [0.1  <  Kn  <  10],  are  characterized  by  the  formation  of  narrow, 

highly  non-equilibrium  zones  (Knudsen  layers)  of  thickness  of  the  order  of  molecular  free  path 
X  ;  the  flow  structure  is  then  determined  by  the  fast  kinetic  processes.  Moreover,  in  case  of 
unsteady  flows,  an  initial  Knudsen  time  interval  is  of  the  order  r0  =  X!  v,  where  v  is  the 

molecular  velocity.  Thus,  the  Knudsen  layer  can  be  computed  accurately  only  by  directly  solving 
the  Boltzmann  equation.  Alternative  approaches  which  approximate  the  Boltzmann  equation  to 
account  for  departure  from  equilibrium,  namely  the  higher-order  continuum  equations  such  as 
the  Burnett  equations  [24,  25],  Grad's  moment  equations  [26]  or  Eu’s  equations  [27]  as  well  the 
particle  methods  such  as  DSMC  [2-4],  have  been  shown  to  have  limitations. 

2.1  Solution  Method  for  the  Classical  Boltzmann  Equation 

In  this  section,  we  briefly  describe  the  finite-difference  method  that  is  currently  employed 
for  solving  the  classical  Boltzmann  equation  (BE)  of  classical  mechanics  for  a  monatomic  gas  in 
translational  non-equilibrium.  In  solving  the  Boltzmann  equation  by  a  finite-diffcrcncc  method, 
the  principal  difficulties  arise  in  calculation  of  the  multi-dimensional  collision  integral;  the 
approximation  to  the  collision  integral  must  tend  to  the  actual  one  as  the  mesh  size  in  the 
velocity  space  tends  to  zero.  In  recent  years,  there  has  been  significant  progress  toward  the 
development  of  an  efficient  and  accurate  numerical  method  for  the  solution  of  Boltzmann 
equation  for  a  monoatomic  gas  [12-15].  In  this  method,  the  Boltzmann  equation  is  solved  on 
fixed  space  and  velocity  grids  by  a  finite-difference  method.  A  projection  method  (that  ensures 
that  the  velocities  before  and  after  collision  belong  to  the  same  grid  of  discrete  ordinates)  is 
employed  for  the  evaluation  of  the  collision  integral  that  ensures  exact  conservation  laws  for 
mass,  momentum,  and  energy  as  well  as  zero  value  of  the  integral  under  thermodynamic 
equilibrium  (when  the  distribution  function  is  Maxwellian).  The  last  property  eliminates  the 
numerical  error  of  computing  the  principal  part  of  the  solution  outside  the  Kundscn  layers  and 
shock  waves  and  thus  considerably  increases  the  accuracy  and  efficiency  of  the  method.  The 
differential  part  of  the  Boltzmann  equation  is  approximated  by  an  explicit  second-order  flux- 
conservative  scheme.  The  combined  system  of  difference  equations  (for  the  collision  integral  and 
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differential  part)  is  solved  by  the  splitting  method  which  splits  the  solution  process  in  two  stages: 
the  collision  relaxation  and  free  molecular  flow'.  The  method  described  below  is  second-order 
accurate  and  is  quite  efficient;  it  is  due  to  Cheremisin  and  his  colleagues  [12-15]. 

Consider  a  gas  of  molecules  without  internal  degrees  of  freedom  and  seek  a  solution  to 
the  2-D  Boltzmann  equation: 


f^+ &  ~+ £v  ~  =  HfJ) s  -L(fJ) + 

ct  ox  oy 

In  equation  (1),  collisions  in  L(f,f)  =  v(f)f  and  (?(/,/)  have  the  form 


(1) 


x  x 

UfJ )  =  j  \  \f-fx-g-b-db-ds-d£x  and  G(f,f)  =  J  J  J /  •  f\  ■  g ■  b  db-de  .  (2) 

-x  0  0  -x  0  0 

A  standard  notation  is  used  in  equation  (2):  b  and  8  arc  the  impact  parameters  of  a 
molecular  collision;  /,/],/  and  fx  arc  functions  of  velocity  vectors  and 

respectively.  The  first  two  vectors  arc  the  pre-collision  velocities;  the  last  two  are  the  post¬ 
collision  ones;  g  =  |£-£,|;  bm  is  the  upper  limit  of  the  impact  parameter;  and  the  vector  £  has 


the  components  and  .  The  space  of  velocity  £  is  restricted  to  a  domain  Q  where  an  N- 


point  uniform  grid  is  defined,  with  grid  points  £,p  and  a  mesh  size  vector /i  =  (/2,,/z-,,/*,) .  Equation 
( 1 )  is  approximated  by  a  set  of  N  equations  for  /„  as: 


%JL+p  H L-r 

dt  ax  ^  dv  ~  p‘ 


(3) 


The  method  used  for  calculating  the  collision  integrals  described  in  [12,  13]  ensures  that 
a  computation  performed  to  any  numerical  accuracy  is  consistent  with  the  conservation  laws.  If 
4 ‘ p  denotes  the  collision  invariant  vector,  4*  =  (!,<*,, ,  then 


IV, -0 

0 

and  at  any  grid  point 


(4) 

(5) 
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Here,  fMp  is  the  equilibrium  Maxwellian  distribution  function  for  gas  molecules  at  a  grid 
pointy .  The  last  condition  substantially  improves  the  accuracy  of  calculation  of  the  integral  /  p 

in  those  subdomains  where  the  state  of  the  gas  is  close  to  the  thermodynamic  equilibrium  [14], 
System  of  equations  (3)  is  solved  by  the  method  of  splitting  with  respect  to  physical 
processes  [15,  16].  On  an  interval [C,C+I],  we  consecutively  solve  the  eollisionless  transport 
equations 


dfp  -  dffi 


s/; 


dt  +^‘ ar+^"‘^r=0’  //=/; 

and  the  eollisional  relaxation  equations 

d/Jdt  =  Ip, 


(6) 


(7) 


written  for  each  grid  point  (the  grid-point  index  is  dropped). 

By  splitting  with  respect  to  spatial  variables,  system  of  equations  (6)  is  approximated  by 
a  second-order  accurate  (in  hx  and  hv )  explicit  seheme  [28].  The  set  of  nonlinear  equation  (6)  is 

solved  by  an  integral  form  [15].  For  time  step  r  =  r'  1  - 


T 


/+! 


fp  +  pdT 

i 


(8) 


To  resolve  fast  kinetic  process,  the  condition  x  <  To  where  To  is  the  mean  time  between 
collisions,  should  be  satisfied.  Cheremisin  and  his  colleagues  have  written  a  Boltzmann  solver 
using  the  above  method.  This  solver  has  been  extensively  validated  by  computing  flows  with 
Knudsen  numbers  ranging  from  0.01  to  10  by  the  author  of  this  report  and  Professor  Cheremisin 
for  computing  hypersonic  flows  with  translational  non-equilibrium  [17].  This  solver  was 
extended  for  solving  the  Generalized  Boltzmann  Equation  (GBE)  for  diatomic  gases;  the 
methodology  is  described  in  Section  2.2. 


2.1.1  The  Stability  Condition  for  the  Algorithm  Described  Above  in  Section  2.1 

The  main  algorithm  of  the  conservative  projection  method  described  above  in  Section  2.1 
has  two  parts.  First,  there  is  the  splitting  of  the  adveetion  equation  and  second  the  integration  of 
relaxation  equation.  For  the  adveetion  equation,  when  an  explicit  1st  order  method  is  applied,  the 
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CFL  condition  defined  by  the  maximum  moleeular  velocity  should  be  satisfied,  i.e. 
^maxrc/ri //;  <  1 .  For  a  higher-order  explicit  method,  instead  of  unity  on  the  right  hand  side  of 

above  condition,  a  number  less  then  1  (for  example  0.5)  appears. 

For  the  solution  at  the  relaxation  stage,  there  are  no  exactly  defined  stability  parameters. 
However  an  important  role  is  played  by  the  parameter £r  that  limits  the  proportion  of  rejected 
contributions.  It  should  be  sufficiently  small.  Usually,  it  should  be  of  the  order  of  E-4  or  E-5.  For 
strict  computation,  it  should  be  equal  to  0,  but  then  the  CPU  time  becomes  very  large.  In 
practice,  the  optimal  choice  of  er  depends  on  a  problem.  The  wrong  choice  may  influence  the 
results  of  the  solution,  and  even  lead  to  the  divergence  of  the  solution.  When  er  is  chosen  and 
the  number  of  integration  (Korobov’s)  nodes  is  given,  the  maximal  value  of  the  time  step  rrat 
the  relaxation  step  can  be  determined.  For  the  stability  of  the  computations,  the  time  step  of  the 
algorithm  should  be  chosen  as  r  =  min(rc/r/,  ,rr) . 

The  Boltzmann  equation  covers  gas  flows  in  all  regimes:  from  free  moleeular,  which  is 
extremely  viscous  from  the  view  point  of  continuum  fluid  dynamics,  to  the  Euler  gas  dynamics 
with  zero  viscosity.  The  algorithm  contains  a  number  of  parameters  whose  values  depend  on  a 
particular  problem.  These  parameters  are:  the  discretization  steps  in  coordinate  and  velocity 
variables,  the  time  step,  and  the  limits  in  the  velocity  spaee.  If  these  parameters  are  not  chosen 
correctly,  the  solution  may  be  wrong  or  diverge,  because  it  doesn’t  correspond  to  the  physical 
reality.  For  example,  if  the  gas  is  very  far  from  the  thermodynamic  equilibrium,  and  one  takes 
the  coordinate  step  h  >  X  (loeal  mean  free  path)  or  the  time  step  is  not  sufficiently  less  the  mean 
moleeular  inter-eollision  time,  the  algorithm  will  produce  wrong  results  and  may  diverge.  The 
use  of  the  same  parameters  for  a  near  continuum  flow  however  may  be  O.K.  The  situation  is 
analogous  to  the  application  of  the  Euler  gas  dynamics  equations  to  a  viscous  low  Reynolds 
flow.  Therefore  the  parameters  of  the  algorithm  that  define  the  approximation  of  a  particular 
problem  should  be  properly  chosen.  The  same  conclusions  are  true  for  the  algorithm  employed 
for  the  solution  of  GBE  or  the  Wang-Chang  Uhlenbeek  equation  described  in  Seetion  2.2 
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2.2  Solution  Method  for  the  Generalized  Boltzmann  Equationn  for  a  Diatomic 
Gas  (Wang-Chang-Uhlenbeck  Equation) 

Nonequilibrium  processes  in  a  gas  with  internal  degrees  of  freedom  of  molecules  can  be 
studied  by  using  the  generalized  Boltzmann  equation  or  the  Wang-Chang-Uhlenbeek  equation, 
where  the  internal  and  translational  degrees  of  freedom  are  considered  in  the  framework  of 
quantum  and  classical  mechanics  respectively.  For  small  deviations  from  local  thermodynamic 
equilibrium,  an  asymptotic  method  developed  for  this  equation  yields  hydrodynamic  equations 
including  the  effect  of  the  internal  degrees  of  freedom  on  the  viseosity  and  thermal  conduction  of 
the  gas  [11].  The  GBE  or  Wang-Chang-Uhlenbeek  equation  (WC-UE)  for  a  diatomic  gas  in 
thermodynamic  non-equilibrium  can  be  written  in  the  form: 


(9) 


In  Equation  (9),  JQ  -  smO df) d(p ,  ft  =  /(/,£,*,/)  is  the  distribution  function,  where  i  is 


the  set  of  quantum  numbers  determining  the  internal  state  of  the  molecule;  is  the  velocity  of 
the  moleeule  in  the  / th  state;  g  =  indices  i,j  and  kj  correspond  to  the  molecular 

states  before  and  after  the  collision  respectively;  and  cr.  is  the  cross  section  for  the  collision 

responsible  for  this  change  of  the  internal  states.  There  is  no  summation  with  respect  to  the 
repeated  index  / . 

The  cross  sections  for  direct  and  inverse  collisions  are  related  as 
ga"(g,0,  (p)d^ld£)j  =g*(Ty(g‘,0,<p)d£kd£l  (10) 

The  magnitude  g'  -|^  -£;|  of  the  velocity  after  the  collision  and  velocities  and 
are  determined  from  the  laws  of  conservation. 

g'=gj =4  -0.5g* ,  and  £,=4  +  0.5g\  (11) 

V  mg 


In  equation  ( 1 1 ),  m  is  the  molecular  mass;  Ae  =  el+ek-el-  e, ,  where  e(.  is  the  energy  of 
the  /th  internal  state  and  ^0-0.5(^i  +  ^l).  The  condition  mg2  >4Ae  determines  the 
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admissibility  of  the  transition (/',  j)->  (k, l) .  We  set  cr*/  =  0  for  forbidden  transitions.  The 
probability  of  the  transition  (i,j)  — >  (kj)  is  defined  as 


CTk! 


Pij  =  —  ,  where  a"  = 
au 

and  satisfies  the  condition  0  <  pkl  <  1  and  the  normalization  condition 

I  rf-1. 


(12) 


(13) 


k.l 

We  assume  that  cr  is  independent  of  the  internal  molecular  state  and  is  equal  to  the 
elastic  scattering  cross  section  cr  =<J0(g,0) .  The  introduction  of  quantities  cr()  and 
obviously  transforms  equation  (9). 


2.2.1  Numerical  Solution  Procedure  for  GBE 


The  problem  of  numerically  solving  the  Wang-Chang-Uhlenbeck  (WC-U)  equation  (29] 
reduces  to  the  construction  of  a  method  of  calculating  its  right-hand  side,  the  generalized 
collision  operator,  which  can  be  represented  as  /,  =  -Lj  +  (J, ,  where 


4  =  cr«  I  j 

j.kJvo 


(14) 


and 

a,  I  JJ MrfgdOdij  ■ 

j,k,iy  n 


(15) 


As  in  the  case  of  the  classical  Boltzmann  equation  (1),  an  effective  numerical  method 
must  ensure  that  the  collision  operator  must  be  (i)  conservative  and  (ii)  be  equal  to  zero  on  the 
equilibrium  distribution  function. 

Operators  I,  and  Gi  are  calculated  on  the  jm  xSn  lattice,  where  Sn  is  the  uniform  lattice 


in  volume  V  of  velocity  space  and  jm  is  the  number  of  quantum  levels. 
Similar  to  the  case  of  a  one-atom  gas,  we  consider  the  functional 

X'  oo 

««>,/)= <t„  Z  fffi Mrfdad&t, 

ij,k,l  J  J  Z. 

-x-xQ 

Taking  functions  0(£)  in  the  fonn 


(16) 
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(17) 


O'  =  0.5 [dnld(^' -4i)  +  dnjd(4’ ~^)\  and 
O*  =  0.5[5„4(4  -  4 )  +  8nld( r  -  £)] ,  (18) 


where  dnq  is  the  Kronecker  delta  and  d(£*— £)  is  the  three-dimensional  delta  function,  we 
obtain 

L„(Z')  =  Q(<t>Ji)  and  (19) 

GXn^QW'Ji)-  (70) 

Functions  (19)  and  (20)  are  calculated  on  the  uniform  cubic  lattice  consisting  of Nv sites 
such  that  (£)v,(aeS0-  Let  (4),.  and  (^), .  .The  values 

(4)v  and  (4)„  for  each  site  of  this  lattice  are  calculated  by  equation  (11).  The  arrangement  of 
vectors  4>  £ . ,  4 ,  and  4  for  thevth  site  of  cubic  lattice  50  is  schematically  shown  in  Fig.  1, 
where  the  three-dimensional  velocity  lattice  is  given  as  a  plane  lattice  and  subscript  v  is  omitted. 
The  value  Ln  v  -  Ln(£y)  calculated  in  site  4  =  g  e  5n  by  equation  (19)  is  determined  by  the  part 


of  the  cubic  sum  for/  =  n ,  av  =  y  and  j\  =  n ,  /?,  =  y  as 

K.y  =  B^8<  8ra,  +  4,Aa  )Av ' 

V 


(21) 


where  B  - 


^cxJ2mV 

N„ 


A,  =  4«. 4a  (a")v^v  sin  4 ,  A,,  =  4fl> 4a (P"),?v  sin  0v ,  and 


fia=f(i,^a,x,t).  In  what  follows,  the  subscript  v  will  be  omitted  where  possible.  Since 
velocities  4  and  4  are  not  in  the  sites  of  lattice  S0,  Gn  is  calculated  with  the  replacement  of 
equation  (18)  by  projector  O"  into  pairs  of  sites  4 ,  gA+s  and  g ,  4  s ,  which  are  nearest  to 
4  and  4  and  are  shown  in  Fig.  13: 


o**(4) = (i  -  r)[dnAd($  - 4) + 44(4  - 4)1 + W'mMt  ~ 4+,) + 4,,-4(4  - 4-,)] .  (22) 

where  s  =  (s,,s2,s3)  is  the  vector  whose  components  take  values  0,-1,  and  1  and  that  determines 
the  site  that  is  nearest  to  4  and  shown  in  Fig.  1.  As  a  result, 

Gn,y  =  BYu  UO  “  K^8nk8yX  +  dnldr,.)  +  r(4*5M«  +  5 n!d y.v-M) >■  •  (22) 
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The  coefficient  rv  is  determined  from  the  laws  of  conserv  ation  for  each  site  of  the  cubic 


formula,  i.e.,  for  each  contribution  A,,  to  operators  Ln  and  GnJ.  The  conservation  of  mass 


follows  from  the  form  of  <P",  and  the  conservation  of  momentum  follows  from  the  symmetric 
arrangement  of  lattice  sites  %x ,  and  ^ with  respect  to  vectors  %k  and  .  In  tenns 


of  the  notation 

+  <  <[+<; 
2  2  2  2 


,  and  £,  = 


2 


(24) 


the  law  of  energy  conservation,  when  contribution  Ay  is  split,  has  the  form  E0  -  (1  -/•)£,  +  rE2 . 


Figure  1:  Schematic  of  an  Inelastic  Collision 

E  -  E 

Therefore  r  =  — 2 - - ,  which  satisfies  the  condition  0  <  r  <  1  because  cither  E,  <  En<  £, 

£;  -  £  102 

or  E2  <  E0  <  £, .  It  is  important  that  rv  is  independent  of  A, .  For  this  reason,  additional  lattice 
sites  ^  ^  +t  ,  and  ^  (  and  coefficient  rv  can  be  preliminarily  determined  for  each 

site  of  cubic  lattice  5, ,  and  then  extended  lattice  S\  can  be  used  repeatedly,  e.g.,  in  various  sites 
of  physical  space.  Each  contribution  Av  can  be  treated  as  the  result  of  a  "collision"  transferring 
molecules  from  sites/,  j  to  sites  A,  //  andA  +  s,  /j-s  .  In  order  to  ensure  that  condition  (ii) 
above  is  satisfied,  we  consider  a  pair  of  inverse  collisions  to  sites  / ,  j  from  sites  A  ,  //  and  A  +  „v, 
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p-s  with  weights  (1-r*)  andr* ,  respectively.  The  probability  p.J  is  determined  from 
equation  (10),  and  coefficient  r  from  the  condition 

£,A(1)(1  -  r)  +  £,  A‘2)  =  £0[AU,(1  -  r)  +  A(  V  ] ,  (25) 


where 

A0)  =  fk,xf,  sin((9)  and  (26) 

A<2>  =  fk.A.+sfi.ti-sPo  S  sin{0) .  (27) 

a"V 

As  a  result,  r  =  — - — - .  (28) 

A"V  +  A,2)(1-a-) 

An  analysis  similar  to  the  above  gives  the  contribution  of  inverse  collisions  in  the  form 

-0(3^  +d„,dm)A {l)  +  r(dnkdrMs  +dnld^  ,)A<2)}V,  and  (29) 

v 

c;,  =  ^A)10-'-')Aro +/--4,!’]}..  (30) 

V 


Finally,  the  collision  operators  are  determined  as 


Kr  =\&JL  +  G,r  =  ^(0^  +  0- 


(31) 


For  any  arbitrarily  lattice  of  integration  sites  Sv,  condition  (ii)  is  satisfied  to  within  an 
accuracy  no  worse  thanO(/?),  where  h  is  the  step  of  lattice^ .  For  degenerate  levels  of  internal 
energy,  equation  (9)  must  be  modified  so  as  to  reduce  the  number  of  substantial  levels. 

Let  jm  levels  be  separated  into  Jm  groups  of  degenerate  levels,  J  =  0,1,...,  JM ,  with 

degeneration  degree  qr  We  determine  the  distribution  function  as F;  =  ^  fa  =  q  .fa , 

where  </6  J  .  Substituting  cr.]  =  <7npfi  into  equation  (9),  summing  this  equation  over  the  groups 
of  levels  i,j, k,  and  /  forming  degenerate  levels  I,J,K  and  L ,  and  returning  to  the  old  notation, 
we  arrive  at  the  equation 


^  +  6  vA  =&  Z  J  jiWjfkfi  -  QkQtf if  j  )gpk,j  j 


(32) 


Oj.kJ-aon 

for  which  equations  (21),  (23),  (29)  and  (30)  are  valid  with  the  change  A  -» qkq,A , 
A0)  ->  qkq, A(,) ,  and  A<2)  ->  qkq, A<2) . 
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2.2.2  Modeling  of  Transition  Probabilities,  Elastic  &  Inelastic  Cross-Sections 

The  computation  of  collision  integral  in  equation  (9)  requires  the  knowledge  of  the 
collision  cross-section  cr*/  responsible  for  the  change  of  internal  states  of  the  molecules.  This 
cross-section  can  be  calculated  from  equation  (12)  if  the  probability  of  transition (/,  y)  —»(£,/) , 
P-1  and  the  clastic  scattering  cross-section  cr(/  are  known  for  a  given  gas  or  a  mixture  of  gases. 

For  DSMC  simulations  of  real  gases,  Koura  [7]  has  given  models  of  clastic  cross-section, 
rotationally  and  vibrationally  inelastic  cross-sections,  and  dissociation  and  recombination  cross- 
sections.  These  cross-sections  can  be  used  for  shock  wave  simulations  in  Nitrogen,  Oxygen,  and 
in  a  mixture  of  Oxygen  and  Nitrogen.  In  the  simulations,  one  may  use  these  models  (or  similar 
models  used  in  DSMC  simulations  [4])  as  well  as  other  models,  which  have  proven  to  be  more 
accurate.  It  has  been  recently  shown  by  Bcylich  [30]  that  Lennard-Jones  (LJ)  model  gives  better 
prediction  for  the  structure  of  normal  shock  waves  in  Helium  and  Argon  than  Hard-Sphere  (HS) 
model  because  the  scattering  behavior  of  LJ  is  very  different  from  the  HS  model.  In  our 
calculations,  we  employ  the  LJ  model  [11]  for  calculating  the  clastic  cross-section  of  diatomic 
gases  e.g.  Nitrogen,  Oxygen,  and  a  mixture  of  Nitrogen  and  Oxygen.  For  calculating  the 
rotationally  inelastic  cross-sections,  we  employ  the  transition  probability  model  due  to  Bcylich 
[22,  23]  based  on  an  interlaced  system  of  rigid  rotors  which  has  shown  excellent  results  for 
shock  wave  structure  in  Nitrogen.  The  model  gives  the  transition  probabilities  as  follows: 

Pf  =  P<M!  K  exp(-A,  -  A2  -  A3  -  A4 ) + — exp(-A3  -  A4 )] , 

«o 

where  A,  =|  Ac,  +  Ae2 1  /  elr0 ,  A2  =  2  |  Ae2  -  Ae,  |  /  eUl, , 

A3  =  4 1  Ae,  |  /(e,r0  +  eri) ,  A4  =  4 1  Ae2  |  /(<?,r0  +  erJ) ,  (33) 

and  <y*'  =  qkq,  !{qiqj ) ,  elr0  =  mg 2  /  4 .  eM  =  elr0  +  en  +  erj . 

In  equation  (33),  qt  is  the  degeneration,  and  e ,  is  the  rotational  energy  of  the  /- th  level.  This 

model  can  be  easily  extended  to  Oxygen.  For  calculating  the  vibrationally  inelastic  cross- 
sections,  we  employ  the  transition  probability  model  described  in  Reference  [7], 
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2.2.3  Computation  of  Shock  Structure  in  a  Diatomic  Gas 

These  computations  were  performed  for  validation  of  the  Generalized  Boltzmann  Solver. 
In  particular,  both  the  computations  with  other  methods  such  as  DSMC  and  BGK-typc  models 
[9]  as  well  as  experimental  data  [31,  32]  are  available  for  shock  structure  in  Nitrogen  for  both 
weak  and  strong  shocks.  The  calculation  of  the  shock  structure  using  equation  (9)  requires  the 
calculation  of  relaxation  of  vibrational  levels  and  rotational  levels.  The  energy  of  j  th  vibrational 

(  l  N\ 


level  is  given  by  e ‘  -tuo 


J  +  ~2 


and  the  equilibrium  distribution  for  the  vibrational  temperature 


Tv  is  given  by 


n]  =  nZv  1  exp 


(34) 


where  Zv  is  the  vibrational  partition  function.  For  Nitrogen, =  3340AT  and  for  Oxygen  it  is 

k 

2230  K,  where  k  is  the  Boltzmann's  constant.  The  energy  level  with  the  rotational  quantum 

ft2j(j+ 1) 


number  j  has  the  degeneracy  of  degree  q  ■  =  2j  +  1  and  energy  e'  =  — 
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,  w  here  Ir  is  the 


moment  of  inertia  of  the  molecule.  The  equilibrium  distribution  of  the  gas  density  n  over  levels 
for  temperature  Tr  is  given  by  the  expression 


n]  =  nZr'(2j  +  l)exp 


KHJ  +  0 

21  T 


(35) 


The  rotational  constant  for  Nitrogen 


2  Ik 


-2.9K  and  for  Oxygen  it  is  2.1  K.  The  values  of 


—j—  and  —  —  for  a  nitrogen  and  oxygen  molecule  have  been  obtained  from  [33].  From  the 

values  of  vibrational  and  rotational  constant  for  nitrogen  and  oxygen,  it  is  obvious  that  the 
vibration  levels  are  excited  at  high  temperatures,  but  the  rotational  levels  are  excited  at  any 
temperature.  Therefore  in  most  of  the  eases,  only  a  few  vibration  levels  can  be  considered.  At 
room  temperature  for  Nitrogen,  about  20  rotational  levels  are  sufficient  for  calculations,  and  at 
least  25  levels  are  needed  for  Oxygen.  When  the  temperature  rises,  the  needed  number  of  levels 


IQ 


increases.  For  shock  wave  in  Nitrogen  at  Mach  2  to  11,  about  30  to  60  levels  may  be  needed. 
Although  the  computational  difficulty  in  number  of  arithmetic  operations  increases  more 
than  J'  times,  fortunately,  for  a  large  number  of  levels  one  can  consider  the  spectrum  as 
continuous  and  use  the  reduced  number  of  efficient  levels.  Such  is  the  approach  applied  in 
statistical  physics  where  the  statistical  sum  is  replaced  by  an  integral.  When  a  reduced  number  of 
levels  are  used,  one  can  consider  the  rotational  levels  in  GBE  or  WC-UE  in  the  framework  of 
classical  mechanics.  In  fact,  Beylich  [22,  23]  has  successfully  employed  this  approach;  he 
computed  the  probability  of  R-T  transfer  for  continuous  spectrum,  and  then  made  it  discrete 
using  the  notion  of  a  rotational  quantum.  One  could  do  the  same,  but  with  an  arbitrary 
“quantum”. 

The  SW  structure  is  formed  as  a  final  stage  of  the  evolution  of  a  discontinuity  in  the  initial 
distribution  function.  The  problem  is  considered  for  the  interval  —  L,  <  x<  L2  with  the 
discontinuity  atx  =  0.  The  initial  distribution  function  on  both  sides  of  the  discontinuity  is 
described  by  the  velocities  and  spectral  levels.  It  has  the  form 


/12(#,x)  =  n'-2[m/(2xT'-2)]m  exp[-  ]^exp(- jfc) , 


2  T 


Qr 


(36) 


where  Qr  denotes  the  statistical  sum.  Parameters  (n,T,u)'2  are  defined  by  the  Rankine- 
Hugoniot  relations  withy  =  7/5.  At  the  boundary,  the  initial  distribution  function  is  kept 
constant.  The  characteristic  dimensional  parameters  needed  in  the  computation  are  the  gas 
density  n,  the  initial  translational  temperature  Tn ,  and  the  mean-free-path  time  r  at  this 


temperature.  The  initial  distnbution  function  at  t  =0  can  be  taken  as  Maxwellian  with 
translational  velocities.  Lattices  Sn  and  Sv  with  about  5000  and  0.5  x  10*1  sites  respectively  are 


sufficient.  Time  integration  is  performed  according  to  the  scheme  in  [12,  14]  using  an 
appropriate  time  step,  e.g.  At  =  0.005r  based  on  the  calculations  performed  in  [17].  In  Section  3, 
the  results  of  computations  for  SW  structure  in  nitrogen  at  high  Mach  numbers  are  presented  for 
R-T  relaxations,  for  both  R-T  and  V-T  relaxations,  and  for  an  inert  mixture  of  two  diatomic 
gases. 
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2.2.4  The  Choice  of  Efficient  Rotational  Levels 

For  real  gases,  except  for  some  light  ones  like  Hydrogen,  the  rotational  energy 
quantum^,  is  very  small  compared  to  the  thermal  energy  (at  the  usual  conditions  when  the  gas 

temperature  is  not  close  to  the  absolute  0  AT).  For  example,  for  Nitrogen  one  has£ro,  Ik  =  2.9 K  , 
and  for  Oxygen  £rot  Ik  =  2.\K  .  Let  the  characteristic  gas  temperature  be T0 .  We  then  define  the 
non-  dimensional  parameter/  =  sml  / kT0 .  When  /  «  1  the  spectrum  contains  a  large  number  of 

levels  and  is  dense.  It  can  be  considered  as  a  “near  continuous”  spectrum.  Because  it  may  be 
difficult  in  some  cases  to  compute  the  real  spectrum  with  all  the  levels,  an  approximate  approach 
can  be  employed  without  loss  of  accuracy.  It  is  not  a  strict  approximation,  but  rather  a 
reasonable  model  of  the  spectrum.  One  considers  the  spectrum  as  continuous  and  applies  some 
step  of  the  discretization  eml  preserv  ing  the  degeneration  rules.  From  physical  consideration,  one 

should  preserve  the  condition  /*  =  ero,  / kTmm  «  1 ,  or  at  least  /*  <  1  .  This  condition  signifies  that 
the  energy  threshold  for  the  excitation  of  the  rotations  everywhere  in  the  flow  remains 
sufficiently  small,  as  it  is  for  the  real  spectrum.  One  may  therefore  suppose  that  better 
approximation  is  obtained  with  small  /  * ,  but  larger  value  of  /  *  leads  to  a  fewer  number  of  the 
efficient  levels  and  may  save  the  computational  effort.  According  to  our  test  computations  of  the 
SW  problem  for  a  wide  range  of  Mach  numbers  from  M=  2  to  25,  the  choice/* <0.25  gives 
quite  satisfactory  results.  Additionally,  it  can  be  noted  from  equation  (33)  that  for  Nitrogen  the 
probabilities  of  transitions  from  the  levels  /,  /  to  the  levels  k,l  don’t  explicitly  contain  the 
value  of  the  rotational  quantum.  This  gives  an  additional  support  for  the  proposed  reduction  in 
the  rotational  levels. 

2.3  Solution  of  Classical  /  Generalized  Boltzmann  Equation  for  a  Mixture  of 
Monoatomic  /  Diatomic  Gases 

For  solving  the  classical  Boltzmann  equation  for  a  mixture  of  monoatomic  gases  or  the 
Generalized  Boltzmann  equation  for  a  mixture  of  diatomic  gases,  these  equations  are  formulated 
in  impulse  space.  These  formulations  are  described  below.  It  should  be  noted  that  these 
formulations  and  solution  methodology'  are  completely  new  and  have  been  developed  for  the  first 
time  by  us  for  3D  Generalized  Boltzmann  equation.  Some  previous  work  for  an  inert  binary 
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mixture  of  monoatomic  gases  for  the  solution  of  ID  classical  Boltzmann  equation  with 
cylindrical  symmetry  in  velocity  space  has  been  reported  by  Raines  [34], 

2.3.1  Mixture  of  Monoatomic  Gases 


The  system  of  Boltzmann  kinetic  equations  for  a  mixture  of  monoatomic  gases  containing 
K  components  is  usually  written  in  the  form 


(37) 


The  collision  integrals  have  the  form 

oo  2  n  bm 

/,=£  J  J  J  (f;f'j-fifj)gbdbd<pd4J,i  =  \,...,KJ  =  \,...,K  (38) 

J  -X  0  0 

Here  bm  is  the  maximum  interaction  distance  and  the  following  abbreviations  have  been  used: 


the  impact  parameters  of  the  binary  collision.  The  six  components  of  the  post  collision  velocities 
vectors  <£  and  are  defined  by  the  three  scalar  conservation  laws  for  the  impulse,  the  energy 

conservation  law,  and  by  the  two  impact  parameters.  For  construction  of  the  conservative 
method  of  evaluation  of  the  collision  integrals  for  a  gas  mixture  one  needs  to  transform  the 
equations  from  velocity  variables  to  the  impulse  variables  defined  asp nr  being  the 
molecular  mass.  Thus, 


(£,x,0->(p,,x,0,  /(£,x, »->/;*( p„x,0. 

From  the  condition  of  normalization  on  the  particle  density  n.  of  a  specie 
tar - j/'rfp  =  ,  one  obtains 


(39) 


/:=»hV, 

In  new  impulse  variables,  the  system  of  Boltzmann  equations  takes  the  form 


(40) 


cl  mj  cx 

The  collision  integrals  take  the  form 


(41) 


7*  =  Z  j  1  J  (fifj  -  f'f'i  )g’bdbd<pdp  J ,  g  =1  p ;.  /  mt  -  p;  /  m,  \ 


(42) 


j  -<c  0  0 


22 


Note  that  in  the  subsequent  equations,  asterisks  in  equations  (41)  and  (42)  will  be  omitted.  The 
integrals  will  be  computed  on  a  limited  Cartesian  impulse  space Q . 

The  system  (41)  is  solved  on  a  uniform  3-dimensional  grid  p  with  N0  points  in  impulse 

spaceQ.  For  brevity,  the  values  of  the  collision  integrals  and  distribution  functions  in  the  grid 
nodes  are  marked  as  L  and  fa  .  The  system  (41 )  of  K  equations  is  transformed  in  the  system  of 

N0K  equations 
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i  =  l . K,y  =  l . N0 


(43) 


dt  >n  dx 

For  evaluation  of  the  integrals  (42),  one  builds  the  8-dimensional  uniform  integration  grid 
s=4i  ,gj  A,<pv  in  impulse  spaces QxQx  Ik xbm  with  Nv  nodes  in  such  a  way  that  impulse 
nodes  pai,  and  pflv  belong  to  the  gridp  .  The  integration  grid  nodes  for  which  the  post-collision 
impulses  p',v  or  pV „ fall  outside  Q  are  excluded. 

Consider  a  mixture  component  /cat  the  grid  nodey.  Introduce  a  combination  of  Dirac  S - 
functions  and  Kroneker  symbols^, ,  with Snl  =\,if  n  =  / ,  and Snl  =0,if  n*l . 


K  =  -  P/.r )  +  4,/(P ,  ~  Pj.r  ) ~  -  P,r )  +  #nAf> 'j  ~  P ,, ) 


(44) 


The  collision  operator  for  the  n  -th  component  at  the  node  yean  be  written  in  the  form 


4r=TZZ  j  j  ]&.?(/,  f  'j-f,  fj  )g  bdbd<pdp.dpj  (45) 

>  /  no  o 

The  integral  is  evaluated  as  a  sum  at  the  gridiE  . 

The  conservative  projection  method  for  evaluation  of  (45)  consists  of  replacing  the  two  last  S  - 
functions  in  (44)  by  decompositions  with  a  splitting  coefficient  rv  <  1  that  has  to  be  defined  from 

the  energy  conservation  law.  For  each  contribution  to  the  integral  sum,  omitting  the  sub-index  v , 
one  makes  the  decomposition 

-  P )  =  Q-r  )<?(PU  -  P,.r )  +  r  <5(pu+,  -  p, ,, ) 

^(p,  — Py.r)  =  C  ~r  —Pj,r)  +  r3(Pj,f,  s  P  j,y  )  (46) 
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In  (46)  the  grid  nodes  pM  and  p  i(are  the  nearest  ones  to  the  post  collision  vectors  p,  and 

p  correspondingly,  whenpu+t,  p  5are  some  complimentary  nearly  located  grid  nodes. 

Hence,  the  contributions  to  the  collision  integral  in  two  near  grid  points  are  replaced  by  the 
weighted  contributions  in  two  pairs  of  the  closest  nodes. 

The  necessary  condition  to  make  this  decomposition  conservative  is  the  fulfillment  of  the 
impulse  conservation  law.  The  transformation  of  the  kinetic  equation  to  the  impulse  space  and 
the  application  of  the  uniform  grid  in  Q  provide  this  condition. 

Let  p,and  p  be  the  impulses  before  collision,  andp'(,  p',.are  the  post  collision  impulses.  We 
select  the  initial  impulses  at  the  grid  nodes,  therefore  they  can  be  presented  as  pj =  k,/r 
andp  =  k2h ,  where  k,  and  k2  are  integer  vectors  and  h  is  the  mesh  of  the  impulse  space. 

The  impulse  conservation  law  gives 


P/  +  P,  =P,  +  P 


(47) 


Let  the  grid  node  nearest  to  p'  be  p  =k,/7  and  that  nearest  to  p'  bep  .  Letpx  =p’-A/7. 
From  (47)  one  gets  pV  =  (k,  +k2  -k,)/7-  Ah ,  therefore  p  =p  +A/7  and  one  gets 


(48) 


P,  +Py  =P^+P/, 


In  a  similar  way,  if  the  nodes  p^+s  and  p^.  are  properly  chosen,  one  can  prove  the  equality 


(49) 


p,  +P,  =P^+,  +P^ 


From  (48)  and  (49)  it  results  that  the  decomposition  (46)  preserves  the  impulse  conservation  law. 
The  decomposition  coefficient  r  can  be  defined  from  the  energy  conservation  law. 

Let  the  v  -th  energy  contribution  to  the  nodes  p(  and  p  be 


(50) 


where  A,,  =  — - f.fiSvK  >  ^  being  the  volume  of  the  Q  space.  Then  the  contribution  to  the 

A  \I  AT  ' 


4^W0 

nodes  p^,,  p  is 
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(51) 


lmt  Zm 

and  the  contribution  to  the  nodes  p^t5 ,  p  is 
E2  =  + 


(52) 


2  mi  2m 

One  should  select  the  nodes  pA+J,  p/win  such  a  way  that  the  value  E0  lies  between  ElandE2. 
The  coefficient  r^can  be  defined  from  the  energy  conservation  law E0  =  (1  - ry)E]  +rx  E2  that 


gives 


r  —  - 


Eq  ~ 

e2-£, 


L,  0  <  r.  <  I 


(53) 


We  have  demonstrated  that  the  transformation  of  the  variables  in  the  system  of  the  Boltzmann 
kinetic  equations  from  velocity  space  to  the  impulse  space  makes  it  possible  to  build  the 
conservative  projection  method  for  the  evaluation  of  the  collision  integrals. 

2.3.2  Mixture  of  Diatomic  Gases 

The  extension  of  the  method  described  in  section  3.3.1  to  a  mixture  of  diatomic  gases  described 
by  a  system  of  Wang  Chang-Uhlenbeck  equations  (WC-UE)  or  generalized  Boltzmann  equations 
(GBE)  can  be  done  in  an  analogous  manner.  The  GBE  for  a  single  component  gas  can  be  written 
in  the  velocity  space  as 

^  r  oc  2x  bm 

dfa 

Px&  i  o'  0 

Here  indices  a,p,x,S mark  the  energy  levels, (oxahp  =  (q^qg)/(qaqfl),  qa  is  the  degeneration  of 
the  energy  level  a ,  is  the  probability  of  transition  from  levels  a,/3  to  the  levels^, 
gap  =|  |.  For  rotational  levels  one  hasga  =2a  +  \  in  most  cases.  For  vibrational  levels,  the 

degeneration  is  absent  andga=l.  We  assume  that  the  degenerations  are  the  same  for  all 
components  of  the  mixture. 

The  generalization  of  (54)  to  the  mixture  of  gases  is  quite  evident: 

Sf,a  e  Of,* 


(54) 


dt 


+£.-^-=L  I J  J 


(55) 


j  faS  -oo  0  0 
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To  build  the  conservative  method  of  evaluation  of  the  collision  operator  on  the  right  hand  side  of 
(55),  the  equation  is  transformed  to  the  impulse  space  in  a  similar  way  as  in  the  case  of  the 
mixture  of  mono  atomic  gases.  One  gets 


of,  a  |  P,  a  df,.a 
dt  m.  dx 


2x  bm 


-£  1 1  I  g«Mbdpdt p,, 


(56) 


fa*  -oo  0  0 

The  important  feature  of  inelastic  collisions  is  that  the  impulse  conservation  law  holds  in  the 
same  form  as  given  in  equation  (47)  for  the  elastic  collisions.  As  a  consequence,  the  similar 
decomposition  of  the  additions  in  the  collision  operator  can  be  made  with  the  single  difference 
that  the  formulae  (50-52)  contain  the  energy  A E  transferred  between  the  translational  and 
internal  degrees  of  freedom.  One  can  exclude  this  value  from  consideration  by  choosing  E0as 
the  kinetic  energy  after  the  collision 


F  —  \  l  P',z  i  P J,s  \ 
Lo  -  A..  •<- —  +  - — ) 
2m  2m 


(57) 


The  formula  (53)  for  the  splitting  coefficient  rv  remains  the  same. 


2.3.3  Special  Case:  Solution  Methodology  for  a  Binary  Mixture  of  Diatomic  Gases 

In  this  section,  we  describe  the  solution  method  for  a  binary  mixture  of  two  diatomic  gases 
as  a  special  case  of  general  methodology  described  in  section  2.3.2.  This  methodology  has  been 
successfully  coded  to  compute  the  hypersonic  shock  structure  in  an  inert  mixture  of  two  gases  as 
shown  in  Figure  1 1.  It  should  be  noted  that  it  is  applicable  to  a  reactive  mixture  as  well.  We 
denote  the  distribution  functions  for  the  mixture  components  as /°(p, x,t),  where  upper  index 

a  marks  the  specie,  and  index/  marks  the  internal  energy  level.  The  generalized  Boltzmann 
kinetic  equations  for  the  two  components  of  the  mixture  can  be  written  as: 


ffi+P_§fl_ 

dt  w,  dx 


=  R"'"  +  R(. 


(1.2) 


¥L+_ pj£ 

dt  m2  dx 


=  R(2-2)  +  R) 


(2.1) 


where  the  collision  operators  are  given  by: 


(58) 
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R 


d.i) 


oo  2itK 

1 1 1  ](/:/;<»“ -f'fX'^gibdbdrt P|J 

y'A7  -oo  o  0 


/?,(2-2)  = 
/?/(,-2)  = 
^:.o  = 


I  j  |  j</,2y>; 

y*/  -oo  o  o 

OO  2xbm 

I I  j 


y*/  -oc  o  o 

oo  2ttb, 


Z  J  \  J  W!/>, 


jkl  _oo  o  0 


-f,2f;)P'2-2>“gilbMd<pllp2J 

-  fi  f,)P'-nt  g.bdbd^j 


These  equations  have  been  written  in  the  impulse  space p  =  {px,p  ,pz)  in  order  to  apply  the 


conservative  projection  method  described  in  section  2.3.2.  The  molecular  masses  should  be 
normalized  by  some  characteristic  mass,  for  examplem0  =2mlm,/(ml  +rn2),  or  simply m0  —  m] . 

For  masses  close  in  values,  the  choice  between  two  methods  of  normalizations  is  not  important. 
We  consider  two  main  cases: 


Casel:  The  characteristic  length  of  the  flow  is  not  large  enough  for  the  vibration  levels  to  be 
excited. 

In  this  case  the  vibration  levels  are  frozen  and  one  may  consider  only  the  energy  transfer 
between  translational  and  rotational  levels  (RT  transfer).  The  collision  operators 
"  include  only  purely  elastic  collisions  and  those  with  RT  energy  transfer. 

For  most  gases,  including  Nitrogen  and  Oxygen,  because  of  small  value  of  the  rotational  energy 
quantum  all  the  collisions  are  not  elastic.  The  characteristic  time  of  the  RT  process  is  4  to  5  times 
larger  than  that  for  the  elastic  relaxation  toward  the  local  thermodynamic  equilibrium.  The  main 
problem  consists  in  evaluation  of  collision  operators R}1^, ft-2’2*, R{0,2\ Rf2,n  .  When  this  is  done, 

the  solution  of  the  system  (58)  can  be  obtained  by  the  application  of  the  usual  splitting  procedure 
described  in  section  3.1  at  a  time  step  r  : 
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(59a) 


Step  (a) 


JC.+.EJC.-0 

dt  m.  dx 


%L+jlW-= 0 

dt  m,  dx 


Step  (b) 


C-//  _  ^(U)  +  ^(1.2) 

dt 


(59b) 


dt 


=  r!2-2)  +  R 


(2.1) 


At  the  second  step  the  collision  operators  are  evaluated  with  the  functions  marked  by  (*).  It  is 
natural  to  apply  at  this  step  a  sub-splitting  in  the  form 

=  *<«..) 


dt 

dt 


=  R 


(1.2) 


(60a) 


(60b) 


and  a  similar  procedure  for  the  second  equation  of  (59b).  Here  in  (60a)  the  collision  operator  is 
computed  by  the  distribution  function  /"  and  in  (60b)  by  the  distribution  function/1”” . 

A  simple  symmetric  modification  wherein  the  equation  (60a)  is  solved  atr/2,  then  (60b)  is 
solved  at  r  ,  and  finally  (60a)  is  solved  at  r  /  2  can  be  applied  for  to  increase  the  accuracy. 

Case2:  The  characteristic  length  of  the  flow  is  so  large  that  the  vibration  levels  are  also  excited. 
In  this  case  one  should  take  into  account  the  huge  difference  in  characteristic  relaxation  times  for 
RT  and  VT  processes.  Let  tvib-sxm^  and  Qvib  denote  the  VT  collision  operator.  Then  the 

inelastic  parts  of  the  corresponding  collision  operators  can  be  estimated  as Qvjh  =sRml,s«  1, 

where  the  indices  on  the  components  and  the  energy  levels  have  been  omitted.  The  kinetic 
equations  then  take  the  form  (vibration  level  is  marked  by  the  index  a  ) 


C  +_pjc 

dt  mx  dx 


-R^'  +  R^+Q-J'  +  Q-f 


(1.2) 


(61a) 
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(61b) 


+  P  Cfi. a  _  ^(2,2) 

dt  m,  5x 


+cJ,+e;.« 


(2,2) 


+  0 


(2.1) 

i,a 


In  principle,  all  the  collision  operators  can  be  evaluated  as  the  GBE  operators,  but  this  task  is 
enormously  difficult  because  the  total  number  of  levels  now  is  equal  to •  amax ,  where  the 

index  indicates  the  numbers  of  the  rotational  and  vibration  levels.  To  make  the  problem  solvable 
for  a  mixture  one  may  apply  the  two-level  model  approach  for  the  RT  collision  operators 
described  below  (without  loss  of  accuracy),  while  the  operators  of  VT  exchange  being  evaluated 
in  the  complete  GBE  form.  In  this  case  the  total  number  of  the  energy  levels  will  be  only  two 
times  the  number  of  vibration  levels  thereby  significantly  reducing  the  computational  effort. 


2.  4  Two  Levels  Kinetic  Model  for  R-T  Relaxation  in  a  Diatomic  Gas 

The  proposed  model  equation  is  aimed  at  simplifying  the  simulation  of  the  rotational- 
translational  (RT)  energy  exchange  in  a  gas.  Such  simplification  is  highly  needed  for  complex 
processes  in  which  rotational  excitation  is  accompanied  by  the  vibration  -  translational  (VT) 
energy  transfer.  The  model  consists  of  2  levels:  the  ground  level  with  the  rotational  energy  e,  =  0 

and  the  excited  level  with  some  energy  e2  >  7max ,  where  Tmay.  is  the  maximum  temperature  in  the 
problem  under  consideration.  We  call  the  proposed  model  as  “2LRT”  model.  The  distribution 
function  is  also  composed  of  two  parts,  /jand /2with  corresponding  populations  of  the  levels 

being  nt  and  n2 .  The  gas  density  is  n  =  w,  +  n2  and  the  rotational  energy  is  Eral  -  e2n2 .  Let  the 
density  of  the  gas  at  some  point  be  77  ,  the  kinetic  energy  Ekm  ,  and  the  rotational  energy  Eml .  One 
can  then  determine  the  populations  of  the  levels  by  the  simple  formulas  n2  -  Ero,  /  e , 
and =77-7?,.  Maximal  value  of  Emt  is  given  by  Eml  -  nTmax ,  therefore  n2  <nTmax  /s2,  and  one 
obtains  0<n2<n  and 77,  >0.  Having Ekm,  one  can  determine  the  equilibrium 
temperature  T  =2(Ekin  +  Erol)/5n,  and  the  equilibrium  rotational 

populations 77,  <  nT .  / e2  ,  77,  =  77  -  772  .  These  parameters  determine  the  equilibrium 

distribution  functions /J  and  f2  M . 
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For  construction  of  the  model  equation  we  begin  with  the  Wang  Chang-Uhlenbeck 
equation  (WC-UE)  for  the  considered  2  levels  system ({i,j,k}  =  1,2) . 

df>ldt  =  Z  j Pu  (fkft  ~  f,  fj  )g,,jbdbd(pd£,J  (62) 

j.k.l 

In  equation  (62),  we  replace  the  collision  operator  by  an  elastic  collision  operator  Qel  and  the 


non-elastic  operator^.  This  replacement  can  not  be  strictly  justified,  as  it  is  assumed  in  [35], 

because  for  RT  exchange  the  purely  elastic  collisions  present  an  exception,  and  about  all  the 
collisions  are  accompanied  with  the  transfer  of  relatively  small  part  of  kinetic  energy  to/from  the 
rotational  energy.  On  the  other  hand,  because  of  small  inelasticity  of  interactions  the  main 
collision  relaxation  process  is  close  to  the  case  of  elastic  collisions,  except  that  one  should  take 
into  account  the  inelastic  transfer  of  the  energy. 

The  elastic  operator  is  the  same  as  the  Boltzmann  collision  integral  for  a  two-component 
gas  mixture: 

Qu,  =  Lf  (f;  ,f'rfifj)g,jbdbd<pdtj  (63) 

j 


The  non-elastic  operator  is  taken  in  a  relaxation  form: 

(64) 

It  was  found  by  a  number  of  numerical  experiments  that  the  choice  for  f’M  in  equation  (64)  as 
the  Maxwellian  distribution  functions  fiM  is  possible,  but  is  not  the  best.  The  function  f'K1 
represents  the  elliptic  distribution  defined  by  the  diagonal  elements  of  the  temperature  tensor 


flu  =  2(0;;0  1  2 ^(-mclUkTl -mc;/2kr  -  me*  /  2kO ,  (65) 


where cx  =  -u,cY  =  £  —v,cz  =£ -w,  and  u,v,w  are  the  components  of  the  bulk  velocity 

vector.  The  components  T‘a  of  the  temperature  tensor  are  defined  by  self-similar  transformation 
of  the  initial  components 

T>TJTq/TkJ  (66) 

The  use  of  the  function  given  in  equation  (65),  instead  of  the  Maxwellian,  means  that  the 
inelastic  operator  preserves  to  some  extent  the  shape  of  the  distribution  function  in  the 
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velocity  space.  The  RT  relaxation  frequency  can  be  defined  as  a  part  of  the  relaxation 
frequency  v  of  the  BGK  model  equation 

vr  =  o,v  (67) 

The  non-elastic  operator  contributes  to  the  evolution  of  the  velocity  distribution  function  toward 
the  equilibrium  state.  To  take  into  account  its  influence  one  should  diminish  the  clastic  collision 
operator  by  a  factor (1  -a2vr),  0<a2  <1 .  Finally,  the  proposed  R-T  relaxation  model  contains 
two  operators,  the  inelastic  operator  given  by  equation  (64)  with  the  frequency  given  by  the 
equation  (67),  and  the  clastic  operator =(1  —  ci1vr)Qii;l.  The  coefficients  gander,  can  be 

determined  from  comparisons  of  the  solutions  of  the  proposed  model  with  solutions  of  the  WC-U 
equation.  This  model  has  been  successfully  applied  to  compute  the  hypersonic  shock  structure  in 
Nitrogen  with  RT  energy  transfer  as  shown  in  Figures  12  and  13.  There  is  good  agreement 
between  the  results  of  two  levels  RT  model  and  those  obtained  with  solution  of  complete  GBE. 
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3.  RESULTS 


Several  milestones  have  been  achieved  which  are  described  below.  The  details  are  given  in 
References  [18-21],  In  what  follows,  some  key  computational  results  are  presented. 

1.  A  3D  Generalized  Boltzmann  Equation  (GBE)  solver  has  been  developed  for  a  Cartesian 
mesh.  The  solver  has  been  validated  by  computing  the  ID  shock  structure  in  nitrogen  for 
Rotational-Translational  (R-T)  relaxations  and  comparing  the  numerical  results  with  the 
experimental  data  for  Mach  numbers  up  to  15.  Figure  2  shows  the  comparison  of  computed 
solution  with  the  experimental  data  of  Alsemeyer  [31]  for  shock  structure  in  nitrogen  at  Mach 
10;  excellent  agreement  is  obtained  [19],  Figures  3,  4  and  5  show  the  details  of  the  computed 
solutions  for  shock  structure  at  M  =15  [19].  Figure  3  shows  the  variation  in  flow  properties  over 
the  shock  thickness.  Figure  4  shows  the  rotational  spectrum  i.e.  the  variation  in  population 
density  of  various  rotational  energy  levels  at  different  locations  in  the  shock  wave,  and  Figure  5 
shows  the  variation  in  rotational  spectral  populations  in  the  shock  region  and  upstream  and 
downstream  of  the  shock.  Figures  4  and  5  arc  very  instructive  in  providing  the  information  as  to 
which  rotational  energy  levels  contribute  most  to  the  shock  structure.  The  solver  has  been 
exercised  successfully  for  computing  the  2D  blunt  body  flows  in  nitrogen  up  to  Mach  7  and  for 
Knudsen  numbers  ranging  from  0.01  to  10  [20].  Figure  6  shows  the  flow  field  contours  for 
density  and  rotational  temperature  at  Kn  =  0.1  when  the  flow  is  in  transitional  regime.  Figures  7 
and  8  show  the  flow  properties  along  the  stagnation  line.  These  results  are  in  good  agreement 
with  the  DSMC  calculations  performed  by  Dr.  Bondar  in  Professor  Ivanov’s  group  in 
Novosibirsk,  Russia  using  the  “SMILE”  code.  The  solver  has  also  been  exercised  successfully 
for  computing  3D  flow  from  a  rectangular  jet  of  nitrogen  exiting  in  vacuum  for  RT  relaxations. 
Figures  9(a)  -  9(e)  show  density,  temperature  and  rotational  temperature  contours  respectively  in 
z  =  0  plane,  and  Figure  9(d)  shows  the  variation  in  density,  temperature  and  rotational 
temperature  along  the  centerline  of  the  3D  rectangular  jet  of  nitrogen  exiting  in  vacuum.  The 
issues  of  stability  of  the  algorithm  and  the  possibility  of  reducing  the  number  of  rotational  levels 
in  the  computations  without  compromising  the  accuracy  of  the  solutions  have  also  been 
rigorously  addressed. 
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Figure  2:  Shock  Wave  Structure  in  Nitrogen  for 
M=10;  n  =  computed  density,  n, experimental 
density,  T  =  total  temp,  Txx  =  translational  temp, 
Trot  =  rotational  temp  (normalized) 
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Figure  4:  Rotational  Spectrum  in  SW  in 
Nitrogen  at  M  =15 
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Figure  3:  Shock  Wave  Structure  in  Nitrogen  for 
M  =15,  n  =  computed  density,  T  =  total  temp, 

Txx  =  Translational  temp,  Trot  =  rotational 
temp  (normalized) 
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Figure  5:  Variation  of  Spectral  Population  in 
SW  in  Nitrogen  at  M=15 
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Figure  6:  Rotational  Temperature  Contours  for  Flow  Past  a  Blunt  Body;  M  =  7 
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Figure  7:  Density  along  the  Stagnation  Line  for  Flow  Past  a  Blunt  Body,  M  =  7 


Kn  =  0.1  Kn  =0.1 

Figure  8:  Rotational  Temperature  along  the  Stagnation  Line  for  Flow  Past  a  Blunt  Body,  M  =7 
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(a)  Density  Contours  in  plane  z  =  0 


(b)  Temperature  Contours  in  plane  z  =0 
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(e)  Rotational  Temperature  Contours  in  z  =  0  (d)  Variation  in  Properties  along  the  Jet  Centerline 
Figure  9:  3D  Supersonic  Flow  of  a  Rectangular  N2  Jet  in  Vacuum,  M  =2,  To  =  200K,  dy  =  2k,  dz  =  8A. 

2.  A  computational  methodology  has  been  developed  to  compute  the  hypersonie  shoek 
structure  in  diatomic  gases  ineluding  both  the  RT  and  Vibrational-Translational  (V-T) 
relaxations  [18].  1-D  shoek  structure  in  nitrogen  has  been  computed  including  both  R-T  and  V-T 
relaxations  and  has  been  validated  by  comparing  the  results  with  experimental  data.  The  V-T 
methodology  is  currently  being  extended  to  3D  in  the  3D  Boltzmann  solver  described  in  the  item 
1  above.  Figure  10  shows  the  variation  in  flow  properties  aeross  the  shoek  at  M  =6  and  M  =10. 


(a)  M  =6  (b)  M  =10 

Figure  10:  SW  in  Nitrogen  including  both  R-T  and  V-T  Relaxations 
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3.  A  computational  methodology  has  been  developed  to  compute  the  hypersonic  shock  structure 
in  a  non-reactive  mixture  of  two  diatomic  gases.  1-D  shock  structure  has  been  computed  in  an  inert 
mixture  of  two  gases  with  R-T  relaxations.  To  accomplish  this,  the  GBE  is  fonnulatcd  and  solved  in 
“impulse  space”  instead  of  velocity  space.  Figure  1 1  shows  the  shock  structure  in  a  binary  mixture 
of  two  gases  of  mass  ratio  ir^/mi  =  2,  and  the  ratio  of  molecular  diameters  d2/di=  1.5.  It  should  be 
noted  that  it  is  easier  to  compute  the  shock  structure  in  an  inert  mixture  of  oxygen  and  nitrogen 
because  the  mass  ratio  is  1 . 143  and  molecular  diameter  ratio  is  1 .0.  The  code  can  be  easily  applied  to 
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Figure  1 1:  SW  in  an  Inert  Mixture  of  Two  Diatomic  Gases 


4.  A  new  two-level  kinetic  model  has  been  developed  for  computing  the  RT  relaxations  in  a 
diatomic  gas  and  has  been  validated  by  comparing  the  results  with  the  solutions  of  complete 
GBE  [1 8]. The  model  is  about  twenty  times  more  efficient  than  the  GBE  in  computing  the  shock 
structure.  It  should  be  noted  that  the  model  is  different  than  the  BGK  model;  it  accounts  for  both 
elastic  and  inelastic  collisions.  Figures  12  and  13  show  the  comparisons  of  flow  properties  using 
the  two-level  kinetic  level  and  the  complete  GBE  for  hypersonic  shock  structure  in  nitrogen  with 
R-T  relaxations  at  M  =  2.4  and  10.  The  agreement  between  the  two  solutions  is  excellent  at 
lower  Mach  number  of  2.4  and  is  reasonable  at  M  =10.  It  should  be  noted  that  a  similar  model  is 
currently  under  development  for  V-T  relaxations.  These  models  can  provide  extremely  efficient 
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solutions  for  gas  flows  of  diatomic  gases  with  both  R-T  and  V-T  relaxations  with  acceptable 
accuracy. 


(a)  Flow  Properties 


(b)  Density 


Figure  12:  Comparison  of  SW  in  Nitrogen  with 
Two-  Level  R-T  Model  (shown  by  *),  M  =2.4 
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Figure  13:  Comparison  of  SW  in  Nitrogen  with  R-T  Relaxations  with  Complete  GBE  and 
Two-Level  R-T  Model  (shown  by  *),  M  =10. 
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